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I 



ABSTRACT 

Stationary solutions of an inviscid and rotational accretion process have been subjected to 
a time-dependent radial perturbation, whose equation includes nonlinearity to any arbitrary 
order. Regardless of the order of nonlinearity, the equation of the perturbation bears a form 
that is remarkably similar to the metric equation of an analogue acoustic black hole. Casting 
the perturbation as a standing wave and maintaining nonlinearity in it up to the second or- 
der, brings out the time-dependence of the perturbation in the form of a Lienard system. A 
dynamical systems analysis of this Lienard system reveals a saddle point in real time, with 
the implication that instabilities will develop in the accreting system when the perturbation 
is extended into the nonlinear regime. The instability of initial sub-critical states may also 
adversely affect the non-perturbative drive of the flow towards a final and stable critical state. 
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1 INTRODUCTION 

Compressi ble fluids, posse ssing angular momentum, will execute rotational motion, when drawn into the gravitational potential well of an 
accretor i Frank et al]|2002h . Devising self-consistent and effective mathematical models to capture the underlying physics of such phenom- 
ena, is a matter of continuing interest to the researcher in astrophysical flows. This is more so true, when the accretor in question is a black 
hole. Astrophysical black holes make their heavy presence felt only by dint of their strong gravitation. Due to their event horizons, no direct 
spectral information about black holes will ever be forthcoming. In this situation, observational evidence for black holes can only be obtained 
indirectly from the way their strong gravitational fields influence any proximate astrophysical fluid. 

To keep matters mathematically simple, a common assumption about rotational accretion is that the fluid flow in question is non-self- 
gravitating (in contrast to, say, a star), and that the motion is axisymm etric. Nevertheless, the m athematical problem remains complex enough, 
involving nonlinear partial differential equations of fluid dynamics dLandau & Lifshital 19871) . One of the principal topics in studies of ax- 
isymmetric accretion is the manner in which angular momentum is transported, enabling a fluid element to move towar ds the accretor in a 
spiral trajectory. Viscosity in the fluid affords the formation of this type of a differentia lly rotating hydrodynamic flow I Lvnden-Belll 19691 : 
Shakura & Sunvaevlll973l ; lLvnden-Bell & Pringlell974l:|Pri ngle 1981 ; Frank et al. 2002), al though the exact nature of the viscous transport is 
a question that has not yet received its final answer ( Prin glell98ll : |Papaloizou & LiiJ 19951 ; Frank et al .l2002h . This, notw ithstanding the effort 
that has already gone into comprehending the role of viscosity in axisymmetric accretion, from multiple perspectives Jshakura & Sunvaevl 
19731: iNovikoy & Thornel| 19731: Evnden-Bell & PringleJ Il974l: lEardlev & LightmarJll975l: Ishapiro et al.l fl97a: Ishakura & Sunvaevl | l976 : 
Ichimarulll977[ Katzll 19771: iBegelmanl 1978l: | Piranlll978l: Liang & Thomsonlll98d:IPringlej 198l | : lBegelman & Meiejl 19821 : Rees et alJll982 : 
Matsumoto et al. 1984; Muchotrzeb-Czernv 1986; Abramowicz et al. 1988; Eggum et al. 1988; Biornsson & Svensson 1991; Luo & Liang 
1994l:lNarayan & Yil 19941 Abramowicz et al.ll995l:lchen et al.ll995l:lchakra barti & T itarchukll995l:lKato et alJl996l:lManmoto et al.ll996l: 
Chakrabarti Hl996iJ:IChen et alJl997l;|Peitz & Applll997l;lBalbus & Hawlevll998l:lFrank et alj2002l;lAfshordi & Paczvnski|2003l:lchakrabarti & Das 



nanian 20( 
alj2009l) . 



2005l :]lslam & Bal busl2005l : IUmurhanetalj200d : (Dasl200llLanzafamel2008l : lsharmal20 08: Subra manian et alj 



2004; Becker & Subramania: 
2008; Bhat tacharjee et 



* soumyab@iitk.ac.in 

f tuhin7048.sg@iiserkol.ac.in 

I arnab.kumar@juet.ac.in 



© 0000 RAS 



2 Bose et al. 



low angular momentum, has been formulated and has become well-established by now (Abramowicz & Zurek 


1981 ; Fukuell987; Chakrabarti 


19891: Nakavama & Fukuell989l: Chakrabartill99C : Kafatos & Yandl994l; Yang & Kafatosl 19951: Parie\Jl996 


Molteni et al J 1996 : Chakrabarti 


1996b; Luetal. 1997; Das 2002; Das et al. 2003; Rav 2003; Barai et al. 2004; Das 2004; Abraham et al. 


2006; Chaudhurv et al. 20061: 


Rav & Bhattachariedl2007al; 


Goswami et alJl2007l:lDas et alj|2007l; iRov & Ravll2009l:lNag et al.ll2012r). In hvdrodvnamic problems it is not 



uncommon to model real fluid flows as inviscid ( Chandrasekh 



lai 



19811) . and as regards rotational accretion, it provides a suitable model for 



describing the flow in the innermost region of the accretion disc, very close to the event horizon of a black hole. Large radial velocities in 
this region imply that the time scale of the dynamic infall process is much smaller than the viscous time scale. Even on large length scales 
the radial velocity has high values due to the fact that the angular momentum of the flow continues to b e low. So when it comes to modelling 
these particular features of axisymmetric accretion, the inv iscid rotational flow certainly has its utility l Igumenshchev & Beloborodovll 19971 : 
iBeloborodov & Illarionov|[200ll ; IProga & Begelmar]|2003l) . 

In feasible models of accretion flows, the boundary conditions require that at large distances from the accretor, the flow is very much 
subso nic, but very close to the accretor, flow ought to b ecome highly supersonic, a fact that is especially true if the accretor is a black 
hole dNovikov & Thorndll973l : lliang & Thomson! 1 1980T) . In other words, in the intermediate region, the bulk flow attains the speed of 
acoustic propagation in the fluid (broadly speaking, the sonic velocity), and becomes transonic in character. Such critical features have 
been investigated in the stationary phase portrait of the flow, and it is known by now that these critical conditions a re the consequences 



of a c oupled first-order dynamical system that can be crafted out of the equations governing the stationary flow (Muchotrzeb-Czern 



1986; iRav & Bhattachariee 



2002 



Afshordi & Paczvriskil liooH; Ichaudhurv et alJbood : IRav & Bhattachariee] l2007al : 



Goswami et alj2007l:lRov & Ray||2007l . l2009l ; lBhattachariee et alj2009l : lNag et alj2012h 



Mandal et al 



zerny 
20071; 



Open critical solutions, connecting the outer boundary of the flow to the surface of an accretor, must pass through saddle points in the 
phase portrait of the stationary solutions, and particularly where the accretor is a black hole, the phase portrait will have mu ltiple saddle points. 



This i s tantamount to saying that an inflow process that is made to traverse all these sa ddle points , must be multi- critical ( Liang & Thomso n! 
1980 | ; | Abramowicz & Zureljl98 JjMuchotrzeb & Paczvnnskill982llMuchotrzebll983l;lMuch^ & Katd 

198S Jchakrabartill 9891,1 199d;lKafatos & Yang 1994l;lYang & Kafatosl 1995l;|Parievll996l : lLu et alJl997l ; |Peitz & Applll99llCaditz & Tsurutal 



as 



& Yang 
sl l2004l ; 



1998; Da: 



2002; Barai et al. 2004; Da: 



Abraham et al. 2006 ; Das et a 1. 2007). a fact that is made possible phys i cally by the occur- 



Nakavama & Fukue 19891; Chakrabartil 199C; Kafatos & Yang 


1994 Yam 


& Kafatos 19951; Lu et alj 1997; Caditz & Tsuruta 1998; 


Das 


2002;; bas et alj|2003l; Ichakrabarti & Das 2004; Abraham et al. 


20061; Das 


2007!: Das et alj 20071: Fukumura & Kazanas! 20071; Lanzafame 


2008; Nagakura & Yamada 2008, 2009). 



Stationary rotational flows may thus be realised in mathematical models by ingeniously tackling the complications that angular mo- 
mentum in the flow presents. For all that, one does not go too far in deriving insight about the dynamic and the nonlinear aspects of the 
fluid flow problem that astrophysical accretion really is. Mathematically, all of fluid dynamics is a nonlinear problem, and accretion flows 
can be no exception to this rule. So it should be a worthwhile exercise to understand the effects of nonlinearity in an accretion process. In 
making a case for nonlinearity, one just needs to look at the well-known nonlinear problem of the realisability of solutions passing through 
saddle points in a stationary phase plot. The very existence of the se critical solutions is threatened b y even an infinitesimal deviation from 
the precisely needed boundary condition to generate the solutions i Rav & Bhattachariedl2002 . 2007a ), and, indeed, the common experience 



regarding the numerical generation of critical solutions in accretion flows is that they are notoriously vulnerable to errors in their outer 
boundary conditions. This difficulty, however, may be overcome by considering the possibility of a non-perturbative temporal evolution of a 
globally sub-critical solution towards the cri tical state, and, in fact, an an alytical study along these very lines has been already been carried 
out, albeit under qualifying approximations (Ray & Bhattacharje e 2007a). Having said this, one needs to realise that nonlinear problems are 
never amenable to an easy mathematical analysis, and, as a matter of fact, there exists no analytical solution of the fully nonlinear coupled 
field equations governing an accretion flow. So in the absence of any analytical formulation of the com plete dynamics of the flow solutions, 
a tentative first step is usually a time-dependent perturbative analysis, and in one such study i Ravlboojl) it was concluded that perturbations 



on the flow do not produce any linear mode with an amplitude that grows in time, i.e. t he background solutions are stable under small 
perturbations. The stability of axisymmetric accretion has been a matter of extensive stud y i Lightman & Eardlevlll974l:lshakura & Sunvaev 



197d:lLivio & Shavivlll977l:IPiranll978l:lKatoll978l:lKato et al.ll98llchen & Taamll993l : lManmoto et al.ll996r . lKato et al.ll996r . lRavl2003l : 



Umurhan & Shavivl l2005; Bhatt achariee & Rav 2007), but conclusions made about stability under a linear regime cannot be projected on 
actual conditions governed by nonlinearity, without courting risks. This work makes an attempt to bridge that gap. 

In this work, a time-dependent, radial perturbation scheme has been implemented and all orders of nonlinearity have been retained in 
the equation of the perturbation that has followed. A most striking feature of the equation of the perturbation is that even on accommodating 
nonlinearity to any arbitrary order, it conforms to the structure of the metric equation of a scalar field in Lorentzian geometry. This fluid 
analogue (an "acoustic black hole"), em ulating many features of a general relativistic black hole, is a matter of pre s ent interest in fluid 
mechanics from diverse po i nts of vi ew (Visser 1998; Schiitzhold & Unruh 2002; Barcelo et a l.l 120051; IVolovik 2005; Si ngha et alj 2005; 
Rav & Bhattachariedl2007<j : lRov & Ravlb007l:lRav & Bhattachariedl2007bl : lDas et alJl2007l : lNag et alJl2012l) . 

The equation of the perturbation is then applied to study the stability of globally sub-critical stationary solutions. Regarding the non- 
perturbative evolution of the accreting system, it is feasible to suggest that the initial condition of the evolution is globally sub-critical, with 
gravity subsequently driving the solution to a critical state, sweeping through an infinitude of intermediate sub-critical states. To accomplish 
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a smooth temporal convergence to a stable critical state, the stability of the sub-critical states is imperative. To investigate this aspect under 
relatively simple nonlinear conditions, all orders of nonlinearity beyond the second order have been truncated in the equation of the pertur- 
bation. Following this, the spatial dependence of the perturbation has been integrated out with the help of well-defined boundary conditions 
on globally sub-critical flows. After this, only t he time-dependent part of the perturb ation is extracted, and, very intriguingly, it acquires the 
mathematical appearance of a Lienard system dStrogat z 1994; Jo rdan & Smithll 19991) . Application of the common analytical tools of dynam- 
ical systems to study the equilibrium features of this Lienard system, shows the existence of a saddle point in real time. This unequivocally 
implies that the stationary background solutions will be unstable in time, if the perturbation is extended into the nonlinear regime. 

To summarise the import of this work, conservative momentum balance and continuity conditions, as appropriate for an inviscid ax- 
isymmetric flow, have been subjected to time-dependent radial perturbations. On allowing nonlinearity to play a part, an instability is seen 
to develop in this hydrodynamic system. And finally, it will not be out of place to me ntion here that th e analytical methods of this work 
are very nearly identical to those of a similar study on spherically symmetric accretion dSen & Ravll2012t) . Considering the great difference 
between the respective geometries of spherically symmetric accretion and axisymmetric accretion, as well as major differences in much of 
the respective physics, the close similarity of the mathematical treatment in both the cases is, indeed, a matter of deep curiosity. The entire 
mathematical treatment so described, and all its attendant physical conclusions, have been presented in what follows. 



2 THE CONDITIONS FOR INVISCID AXISYMMETRIC ACCRETION 

In models of accretion discs, imposing the condition of hydrostatic equilibrium along the vertical direction JMatsumoto et al J 1984 Fran k et al 
l2002h and also performing a vertical integration, will effectively collapse the vertical geometry of the flow on the equatorial plane of the disc. 
The equatorial flow is described by two coupled fields, the radial drift velocity, v, and the surface density, E, of which, the latter is defined 
by vertical ly integrating the volume density, p, over the disc thickness, H. This gives E = pH, and in terms of E, the continuity equation 
will go as jFrank et al.ll2 002) 



9s + ^(H = o. 



(1) 



dt r dr 

The axisymmetric accretion flow is driven by the gravitational field of a centrally located non-rotating black hole, but the structure of 
the neighbouring geometry, through which the flow takes pl ace, is still set according to the Newtonian construct of space and time, by taking 
recourse to what is known as a pseudo-Newtonian potential l lPaczvnski & W iita 1980; Nowak & Wago ner 1991 ; Artemova et al.l 19 96"). This 
effectively substitutes the properties of the Schwarzschild space-time geometry with a potential function, and to a large extent the analytical 
results pertaining to the rotational flow remain unaffected by the choice of a particular pseudo-Newtonian potential. Involving such a potential, 
the height of the disc, under hydrostatic equilibrium in the vertical direction, is given as H = r y~ 1 ^ 2 r(a/vK), with the local speed of 
sound, a, and the local Keplerian velocity, vk, defined, respectively, by a 2 = 7-P/p and «| = r$' (here the prime indicates a derivative with 
respect to r). The pressure, P, as it has been introduced in the definition of a, is expressed in terms of the volume density, p, by a polytropic 
equation of state, P = kp~' , with 7 being the polytropic exponent. In consequence of this definition of P, the speed of sou nd may also be 



noted from a 2 = dP/dp = -ykp 1 1 . It can be shown mathematically, by going back to the first law of thermodynamics dChandrasekhai 



1939), that 7 varies between unity (the isothermal limit) and cp /cv , which is the ratio of the two coefficients of specific heat capacity of a 
gas (corresponding simply to the conserved adiabatic limit), i.e. 1 ^ 7 ^ cp/cy. So the polytropic prescription is of a much more general 
scope than the simple conserved adiabatic case, and is suited well for the study of open systems like astrophysical flows. 

Now, using the relationship between a and p, the disc height can be explicitly written in terms of the standard fluid flow variables as 



H = ( 7 fc) 



1/2 pCr-D/V/* 

vw 1 



and with the foregoing result, the continuity condition, as given by equation Q}, can be recast as 



dt V \ 



+ 



*3/2 flj. 



p (7+l)/2 w 3/2 



= 0, 



(2) 



(3) 



which is one of the two required mathematical conditions to govern the dynamics of the coupled fields, v and p, in the radial direction. 

To ascertain the second condition necessary for determining the dynamics in the radial direction, one first has to look a t the dynamics 
along the azimut hal coordinate. This gives the condition for the balance of the specific angular momentum of the flow as 1 Pringle|[l98ll ; 
Frank et alj|2002h 



Qt ' r dr 1 J » A r 



2nr \ dr 



in which, D. is the local angular velocity of the flow, and the torque, Q, is to be read as 
g = 2^r 2 f^j , 



(4) 



(5) 
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wi th v being the kinematic v iscosity associated with the flow. Using the continuity condition, as equation (TJ gives it, and g oing by 



thelshakura & Sunvaev 1 1973h prescription for the kinematic viscosity, v = aaH, it would be easy to reduce equation (0 to the form jFrank et"al 
200ilNaravan & Yill994h , 

on 



id_ 

V dt 



<:) 



1 d 
pvrH dr 



apHa 2 r 3 



dr 



(6) 



with S1k being defined from vk = tQk- In equation (6J, the condition for an inviscid and axially symmetric rotation of the flow is achieved 
by putting a = 0. This delivers the simplest possible integral solution, r 2 Q, = A, with the constant of the motion, A, being interpreted 
physically as the globally constant specific angular momentum of the flow. This interpretation affords great convenience in setting down the 
second mathematical condition of the flow along the radial direction. This is the condition of the radial momentum balance (in essence the 
conservative Euler equation), in which the centrifugal term, arising due to the rotational motion of the flow, is fixed as A 2 /r 3 . After that, the 
radial momentum balance equation is written as 



— +v— + -— + $'(r) - — - 
dt dr p dr r s 



(7) 



On specifying the functions, $(r) and P, equations l[3j and 10 can give a complete hydrodynamic description of the axisymmetric 
flow in terms of the two fields, v(r, t) and p(r,t). From these dynamic variables, the steady solutions of the flow are obtained by making 
explicit time-dependence disappear, i.e. d/dt = 0. The resulting differential equat ions of the inviscid rotational flow, involving full spatial 
derivatives only, can then be easily integrated to get the stationary global solutions dChakrabarti|[l989l ; lDasll2002ljDas et alj|2003h . A notice- 
able feature of these stationary solutions is that they remain invariant unde r the transformation v — > — v, i.e. the mathematical problem of 
inflows (v < 0) and outflows (v > 0) is identical in the stationary state dRav & Bhattach ariee 2007a). This invariance has some adverse 
implications for the critical flows in inviscid axisymmetric accretion processes. Critical solutions pass throu gh saddle points in the stationary 
phase portrait of the flow l Ray & Bhattacharieej [2002; Cha udhurv et al.|[2006l : Ray & B hattachar ieell2007al) . and it can be demonstrated that 
generating a stationary solution through a saddl e point will be impossible by any physical means, because it calls for an infinite precision 
in the requ ired outer boundary cond ition dRav & Bhattacha riee 2002, 2007a). Nevertheless, criticality is not a matter of doubt in accretion 
processes lliang & Thomsonlll98Cl) . The key to resolving this paradox lies in considering explicit time-dependence in the flow, because of 
which, as one may note from equations $3$ and Q, the invariance under the transformation, v — > —v, breaks down. Obviously then, a 
choice of inflows (v < 0) or outflows (v > 0) has to be made at the very beginning (at t = 0, as it were), and solutions generated thereafter 
will be free of all the difficulties associated with the presence of a saddle point in the stationary flow. 

In the stationary phase portrait of the flow, criticality is obtained when the bulk flow velocity matches the speed with which an acoustic 
wave can propagate through the moving compressible fluid. In the case of an inviscid, polytropic, rotational flow, this speed is not ex actly 
given by the sonic condition, but differs slightly from it by a constant numerical factor, \[2 (7 + 1) -1 ^ 2 , due to the geometry of the flow 1 Ray 
l2003l:lchaudhurv et al.ll2006l : iRav & Bhattacharjeell2007al : lRov & Ravll2009l : lNag et al.ll2012h . For the conserved, inviscid flow, the critical 
points are either saddle point s, through which open solutions may pass, or they are centre-type points, around which the stationary solutions 
will form closed trajectories Jchaudhurv et alfeood) . Dependin g on the global nature of the phase portra it, the solutions which pass through 
the saddle points can be either homoclinic or heteroclinic Jchakrabartilll989l ; lDasll20o"2l ; lDas et al JI2OO3I) - The number of critical points that 
may be obtained, depends on the choice of the potential, 4>(r). For the simple Newtonian potential, only two critical points result dRav 



2003), one a centre-type point and the other a saddle point. On the other hand, in studies of axisymmetric accretion on to a black hole, it is an 



expedient practice to employ a pseud o-Newtonian potential to drive the flow. In that case, the number of critical points will exceed two, and 
going by the nature of critical points 1 Jordan & Smith|[T999l) , it happens that there will be multiple saddle points. This can have a significant 

bearing on the way a fluid element may reach the accretor, after having started from an outer boundary point. 

I mposing various boundary conditions on the stationary integral solutions, will result in multiple classes of flow Jchakrabartil 19891 : IPasl 
200i iDas et al.ll2003l) . Of these, the one of enduring interest in rotational accretion obeys the boundary conditions, v — > as r — > 00 



(the outer boundary condition) and v > a at small values of r. The inner boundary condition naturally suggests itse lf when it comes to 
accre t ion onto a black hole, fo r which, the final infall across the event horizon must necessarily be highly supersonic dNovikov & Thorne 
19731: lliang & Thomsorj[l98(j|) . Therefore, it is quite obvious that an open solution that starts with a very low velocity far away from the 



accretor, but has to allow a fluid element to reach the accretor with supersonic speeds, perforce has to pass through a saddle point (where 
criticality in the flow is attained). This implies that travelling from an outer boundary point, the first critical point that a physically relevant 
open inflow solution encounters is a saddle point. If, however, between this point and the surface of the accretor, there are other saddle points 
to be encountered, then what transpires is a multi-critical flow. In such a flow, a fluid element will reach the accretor aft er having travelled 
through more than one sad dle point, and in between two successive saddle points, the flow is bound to suffer a shock Jchakrabartill 19891 : 
Dasll2002l ; IDas et al ] |2003l) . It is at the discontinuity of a shock, that a solution is demoted from its super-criti cal state to a sub-critical one , 



following whic h, the solution again has to regain super-criticality by travelling through another saddle point Jchakrabarti|[l989l: iDaslbooi 
Das et al]|2003l) . Evidently, multi-criticality and shocks are closely related to each other in models of conserved, rotational inflows onto a 
black hole. 

Thus far, the properties of the accretion flow have been understood from a stationary perspective, which is relatively easy to follow. It 
is the dynamics of the flow (even a conservative one) that offers a pr oblem of greater mathematical complexity. In comparatively simpler 
studies involving time-dependence (Ray 2 0031 ; Chaudhurv et al. 2006), the inviscid and axisymmetric flow is found to be stable under the 
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effect of linearised perturbations. However, this does not shed much light on the way the velocity and density fields in an accretion process 
may evolve non-perturbatively in time. In such a mathematical problem, one is required to work with a coupled set of nonlinear partial 
differential equations, as implied by equations I0 and ifTJ, but no analytical solution exists for these coupled, dynamic nonlinear equations. 

Now, so far as generating the critical flow is concerned, the non-perturbative dynamic evolution of global v(r, t) and p(r, t) profiles is 
very crucial. It is the way in which the two fields evolve vis-a-vis each other that decides if the critical state would be achieved or not. The 
dynamic process should be envisaged mathematically as one in which initially both the radial drift velocity and the density fields, v(r, t) and 
p(r, t), are sub-critical and nearly uniform for all values of r, in the absence of any driving force. Then with the introduction of a gravitational 
field (at t — 0), about whose centre, the fluid distribution may randomly possess some net angular momentum, the hydrodynamic fields, v 
and p, start evolving in time. In the regions where the temporal growth of v outpaces the temporal growth of p (to which a is connected), and 
gravity (going as r -2 ) dominates the inhibitive centrifugal effects of rotation (going as r~ 3 ), the infall process will become super-critical. 



Other wise, it will continue to remain sub-critical. Under the approximation of a "pressureless" motion of a fluid in a g ravitational field ( Shu 



1991), qualified support for the attainment of criticality was provided from a non-perturbative dynamic perspective (Ray & Bhattachariee 



2007a), guided by the physical criterion that the total specific mechanical energy at the end of the evolution will re main the same as what it was 



at the start o f the evo lution. A si milar principle makes transonicity possible in spherically symmetric accretion (Bondi 1952; Garlick 1979; 
Ray & Bhattacharieel2002l ; lRov & Ravh007l) . Nevertheless, a completely analytical solution of the dynamic nonlinear problem continues to 



be elusive. 



3 NONLINEARITY IN THE PERTURBATIVE ANALYSIS 



Equations l[3) and are easy to integrate in their stationary limits, and the resulting velocity and density fields, derived from these two 
equations, have only spatial profiles, v = vq (r) and p = po(r). The meth od of the perturbative analysis that is adopted here, and imple- 
mented originally in spherically symmetric accretion IPetterson et al. 1980h . entails applying small time-dependent radial perturbations on 
the stationary profiles, vo(r) and po(r), and then linearising the perturbed quantities. This, however, does not offer much insight into the 
time-dependent evolutionary aspects of the hydrodynamic flow. The next logical act, therefore, is to incorporate nonlinearity in the perturba- 
tive method. With the inclusion of nonlinearity in progressively higher orders, the perturbative analysis incrementally approaches the actual 
time-dependent evolution of the global solutions, after it has started with a given stationary profile at t = (it makes physical sense to 
suggest that this initial profile is spatially sub-critical over the entire flow domain). 

The prescription for the perturbation is v(r, t) = Wo(r) + v'(r, t) and p(r, t) = po{r) + p'(r, t), in which the primed quantities indicate 
a perturbation about a stationary background. It is now necessary to define a new variable, f(r,t) = p^ +1)/2 vr 3/2 /y/W, following a 



similar mathematical proced ure employed in some previous studies on inviscid, axisymmetric accretion (Rav 20q3; Chaudhurv et al . 2006; 
Rav & Bhattacharieal2007a) . It is easy to see that this variable emerges as a constant of the motion from the stationary limit of equation ((5). 
This constant, /o, can be identified closely with the matter flow rate, within a fixed geometrical factor, and in terms of Vo and po, this constant 
is given as fo = p 1+1 ^ 2 vor 3 ^ 2 /\/W' '. On applying the perturbation scheme for v and p, the perturbation in /, without losing anything of 
nonlinearity, is derived as 

f ft 2 po Vq P 2 po V ' 

in which, 



C = l + 



1 • 2 



P_ 

Po 



1-2-3 



+ • 



(8) 



(9) 



and f3 = \/2(7 + 1) 1//2 - Equation <[Sj connects the perturbed quantities, v' , p' and /', to one another. To get a relation between only p' and 
/', one has to go back to equation 10, and apply the perturbation scheme on it. This will result in 

(10) 



(11) 



i. 






d 


( f'\ 


dt 


\P 2 


po) 




\fo) 



To obtain a similar relationship solely between v' and /', one needs to combine the conditions given in equations ([8) and l |10t , to get 

dv' 



df df 
dt dr 



_ v 
~dt~7 

It is worth stressing at this point that in equations {8]l, d 10b and (TTJ, all orders of nonlinearity have been maintained. Adhering to the same 
principle, perturbing equation Q yields 

dv' 



h _9 

dt dr 



vov + 



+ 



al 



7' 



i5> 

3=1 



P P0 



0. 



in which ao = ao(r), is the steady value of the local speed of sound, and 



7-1 
.7+1 



7+1 



7+1 



7+1 



(J - 1) 



(12) 



(13) 
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Now taking the second partial time derivative of equation d!2t . and making use of equations JlOb . i ll It and the second partial time derivative 
of equation i ll It . one obtains a fully nonlinear equation of the perturbation, running in a symmetric form as 

£ U u %) + 4 (h»¥) + 1 (h"°£) + 4- (irg.) = o, d4) 



dt \ dt J dt \ dr J dr \ dt J dr \ dr 
in which, 

Going by the symmetry of equation J 1 4b . it can be recast in a compact form as 

9m (/»""&/') = 0, (16) 

with the Greek indices running from to 1, under the equivalence that stands for t, and 1 stands for r. At this stage it should be important to 
realise that equation J16b , or equivalently, equation J14b . is a nonlinear equation containing arbitrary orders of nonlinearity in the perturbative 
expansion. All of the nonlinearity is carried i n the metric elements, W v . If one were to have wo rked with a linearised equation only, then 
/i M " could be read from the symmetric matrix l Chaudhurv et al.| [2006: Ray & Bhattacharjee 2007^), 



h M = -r 2 2 2 • ( 17 > 



_ wo / 1 wo 

Jo V ^0 «o - P «0 

Now, in Lorentzian geometry the d'Alembertian for a scalar field in curved space is expressed in terms of the metric, g M „, as 

= ^==d M (V=i gTduv) , (18) 



with being the inverse of the matrix implied by <7 M „ JvisseJl998ljBarcel6 et al.l2005l) . Comparing equations l !16b and l !18b with each other, 



one could look for an equivalence between W and sf—gg 1 *". What can easily be appreciated from this comparison is that equation d 1 6b 
gives an expression for /' which is of the type given by equation d 1 8b . In the linear order, the metrical part of equation l !16b . as equation l !17b 
shows it, may then be extracted, and its inverse will incorporate the notion of the horizon of an acoustic black hole, when vq — f3 2 a%. 
This point of view has features that are similar to the metric of a wave equation obtained by setting the velocity of an irrotatio nal, inviscid 
and barotropic flui d flow as the gradient of a scalar potential, and then by imposing a perturbation on this scalar potential I Visser 19981 ; 



Barcelo et alj|2005l) . In contrast to this approach of exploiting the conservative nature of the flow to craft a scalar potential, the derivation of 
equation d 1 6b makes use of the continuity condition. The latter method is more robust because the continuity condition is based on matter 
conservation, which is a firmer conservation principle than that of energy conservation, on which the scalar-potential approach is founded. If, 
for instance, the fluid system under study were to have contained dissipating factors (as it is happens in models of axisymmetric accretion), 
then the potential-flow formalism would have been ineffective, but the present method of making use of the matter conservation principle, 
would still have delivered an equation of the perturbation. 

However, the close correspondence between the physics of acoustic flows and many features of black hole physics is valid only as far as 
the linear ordering goes. When nonlinearity is to be accounted for, then instead of equation d 1 7b . it will be equations d 1 5b which will define 
the elements, W v , depending on the order of nonlinearity that one wishes to retain (in principle one could go up to any arbitrary order). The 
first serious consequence of including nonlinearity is to lose the argument in favour of an inflow solution attaining the critical status, because 
the descripti on of b? v , as stated in equation < !17b . will not suffice any longer. This view is in perfect conformity with a numerical study 
conducted bv llVlach & Maled J2008I) for the case of spherically symmetric accretion, in which it was shown that if the perturbations were to 
become strong then the analogy between the "sonic horizon" and the event horizon of a black hole would not hold. More to this point, when 
discussing an axisymmetric flow it has been pos sible to invoke th e critical aspects of spherical geometry, because the mathematical methods 
employed in both the cases are nearly identical ISen & Ravllioii) . While this unifying principle is worthy of a serious note in its own right, 



another very remarkable fact to have emerged in consequence of including nonlinearity in the perturbative analysis is that regardless of the 
order of nonlinearity that one may desire to go up to, the symmetric form of the Lorentzian metric equation will remain unchanged, as shown 
very clearly by equation d 1 6b - For t he laboratory fluid problem of the hydraulic jump, a similar type of symmetry is known to exist, going up 
to the second order of nonlinearity jRav & Bhattac harjee 2007tJ). 



4 STANDING WAVES ON SUB-CRITICAL INFLOWS 



All physically relevant stationary inflow solutions obey the outer boundary condition, v(r) — > as r — > 00. Among these solutions, a 
critical inflow will reach the accretor with a high super-critical speed, but somewhere along the way, this flow will also undergo a discontinuity 
due to a shock. So critical inflows are highly su b-critical at the outer boundary, highly super -critical near the accretor (the inner boundary), 
and have a discontinuity in the interim region 1 Chakrabarti 19891 : Dasll2002l : Das et al.ll2003f) . On the other hand, there is an entire class of 
inflow solutions which are globally sub-critical, conforming to the inner boundary condition, v(r) — > as r — > 0. From the point of view 
of a gravity-driven evolution of an inflow solution to a critical state, the sub-critical flows have a great importance, because the initial state 
of an evolution, as well as the intermediate states in the progression towards criticality, should realistically be sub-critical. So the stability 
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of globally sub-critical solutions will have a significant bearing on how a critical solution will develop even tually. Imposing an Eulerian 



perturbation on these sub-critical inflows, their stability has been studied by now under a linearised regime (Ray 2003; Chaudhury et al 



2006), and the amplitude of the perturbation, fashioned to be a standing wave in that case, was seen to maintain a constant profile in time. 



Further, the perturbation was also made to behave like a travelling wave under a WKB approximation, and it was found that although the 
relative amplitude of the perturbation grew locally when the bulk velocity decreased, the total energy in the wave remained fixed. In all of 
these respects one may then say that the solutions do not exhibit any obvious instability. However, it is never prudent to extend this argument 
too far, especially when one considers nonlinearity in the perturbative effects, as it rightly ought to be done in a fluid flow problem. 

Now a nonlinear equation of the perturbation, that accommodates nonlinearity up to any desired order, has already been derived, as 
equation dl4b shows. This equation can then be applied to study the stability of stationary sub-critical flows in a nonlinear regime. The 
perturbation is designed to behave like a standing wave about a globally sub-critical stationary solution, obeying the boundary condition 
that the spatial part of the perturbation vanishes at two radial points in the axisymmetric geometry — one at a great distance from the 
accretor (the outer boundary), and the other very close to it (the inner boundary). While the former boundary is a self-evident fact of the 
flow, there is a certain measure of difficulty in identifying the latter. The guiding principle behind the choice of the two boundaries is that 
the background stationary solution should be continuous in the interim region. For a globally continuous sub-critical solution, that connects 
the outer boundary to the accretor, the inner boundary is obviously the surface of the accretor itself. If, on the other hand, even a sub-critical 
inflow is also disrupted by a shock, then the inner boundary should be the standing front of the shock itself. It is conceivable that no part 
of the perturbation on the background flow may percolate through the shock front, and so, the discontinuous front itself may be set as a 
boundary for t he perturbation. Such piecew ise continuity of a stability analysis, on either side of a discontinuity, is not uncommon in studies 
on fluid flows jRav & Bhattacharieell2007bl) . 



The mathematical treatment involving nonlinearity will be confined to the second order only (the lowest order of nonlinearity). Even 
simplified so, the entire procedure will still carry much of the complications associated with a nonlinear problem. The restriction of not going 
beyond the second order of nonlinearity implies that h ltu in equations j 1 5b will contain primed quantities in their first power only. Taken 
together with equation d 14b . this will preserve all the terms which are nonlinear in the second order. So, performing the necessary expansion 
of v = Vo + v' , p — po + p and / = fo + f in equations J 1 5b up to the first order only, and defining a new set of metric elements, 
(f v = one obtains 

{q^d,j') = 0, (19) 

in which p and v are to be readjust as in equation ( 1 161 ). In the preceding expression, the elements, g^ 1 ", carry all the three perturbed quantities, 
p' , v' and /'. The next process to perform is to substitute both p' and v' in terms of /', since equation dl9b is over /' only. To make this 
substitution possible, first one has to make use of equation ([8} to close v' in terms of p' and /' in all g M ". While doing so, the product term of 
p' and v' in equation ([8j is to be ignored, because including it will raise equation d!9b to the third order of nonlinearity. By the same token, 
one also has to take only £ = 1 from equation {9}. Once v' has been eliminated in this manner, one has to write p' in terms of /'. This can be 
done by invoking equation jj 1 0b . with the reasoning that if p' and /' are both separable functions of space and time, with the time part being 
oscillatory (all of which are standard mathematical prescriptions in perturbative analysis), then 

4-=^«T> (20) 

P PO JO 

with a being a function of r only (which extends a crucial advantage in simplifying much of the calculations to follow). The exact functional 
form of a(r) will be determined by the way the spatial part of /' is set up. It is known that a(r) is indeed a real functio n, dependent on vo 
and ao, when the spatial part of /' is set down as a power series in the WKB approximation (Rav 2003; Chaudhury et al . 2006). 
Following all of these algebraic details, the elements, g ,u/ , in equation d 1 9b . can finally be expressed entirely in terms of /' as 

tt ( , j-tt f | tr 2 f , r tr f \ rt 2 f , /-rt f \ rr I 2 a 2 2\ , 3 /-rr f ,~ , , 

q =vo\l + e£ = v o I 1 + e C t) > 9 = v o I 1 + ~r J > 1 = v (v - /3 a ) + ev £ — , (21) 

in all of which, e has been introduced as a nonlinear "switch" parameter to keep track of all the nonlinear terms. When e = 0, only linearity 
remains. In fact, in this limit one converges to the familiar linear result implied by equation l !17b . In the opposite extreme, when e = 1, in 
addition to the linear effects, the lowest order of nonlinearity (the second order) becomes activated in equation l |19b . It may be noted that 
equations J2 1 b also contain the factors, , all of which are to be read as 



/-tt /-tr /-rt -< n /- 



3 + 



7 -3\ p 2 al 



(22) 



7+1, 

with £' T having been determined by going up to j = 2 in equation d 1 3b . Taking equations l |19b , ( 121b and ( 122b together, one gets a nonlinear 
equation of the perturbation, completed up to the second order, without the loss of any relevant term. 

To render equation d 1 9b . along with all g M " and £^ v , into a workable form, it will first have to be written explicitly, and then divided 
throughout by vq. While doing so, the symmetry lent by £' r = £ is also to be exploited. The desirable form of the equation of the 
perturbation should be such that its leading term would be a second-order partial time derivative of /', with unity as its coefficient. To arrive 
at this form, an intermediate step will involve a division by I + e£ if' /fo), which, binomially, is the equivalent of a multiplication by 
1 — e£ (f /fo), with a truncation applied thereafter. This is dictated by the simple principle that to keep only the second-order nonlinear 
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terms, it will suffice to retain just those terms which carry e in its first power. The result of this entire exercise is 



o 2 f , _ d ( Of 



1 d 

+ —IT 

vo or 



i 2 
^0 [v 



a 2 2\ 9f 

' c?r 



e JtttlWY ,3 ( t df' 2 



2 dr dt 



2vo dr \ dr 



0f tt fl d ( Of 

25 ; Tr r ST 



vo dr 



I 2 
vo (v 



r 2 2 \ d f 
■ P a o) -5- 
' dr 



(23) 



in which, if one were to set e = 0, then what would remain would be the linear solution discussed in detail in some earlier works I Ravi 20031 : 
Chaudhurv et"al]|2006h . To progress further, a solution of f' (r, t), separable in space and time, is to be applied. This will bear the form, 
f'(r,t) — R(r)<j>(t). Using this separable solution in equation d23t . then multiplying the resulting expression throughout by voR, and then 
performing some algebraic simplifications by partial integrations, will finally lead to 



4>v R 2 + <j>4- (voR) 2 +<t>\4- 
ar dr 



vo I 2 a 2 2\ d-R 

7 O*>-0°o)-dT 



, 1 :l 1 , I d-R 
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+ T 
Jo 



b 2 £ u voR 3 + < 



voR)+Z y^r-« R ^(^ R ) 
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dr 



+ <t> 2 1 vo {yl 



p 2 a 2 ) — — (f 
dr dr 
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■ 2 
vo 



a 2 2- 

- P a 



dR* 
dr 



+ 



dr 



. v% dR 3 
3 dr 



0,(24) 



in which the overdots indicate full derivatives in time. Quite evidently, equation d24t is a second-order nonlinear differential equation in both 
space and time. The way forward now is to integrate all spatial dependence out of equation d24t . and then study the nonlinear features of the 
time-dependent part. The integration over the spatial part will necessitate invoking two boundary conditions. One of these (the outer boundary 
condition) is situated very far from the accretor, while the other one (the inner boundary condition) could either be very close to the accretor, 
or at a standing shock front where the background solution becomes discontinuous. At both of these boundary points, the perturbation will 
have a vanishing amplitude in time, while the background solution will maintain a continuity in the interim region. The boundary conditions 
will ensure that all the "surface" terms of the integrals in equation d24t will vanish (which explains the tedious mathematical exercise to 
extract several "surface" terms). So after carrying out the required integration on equation J24b . over the entire region trapped between the 
two specified boundaries, all that will remain is the purely time-dependent part, having the form, 



+ e[A4> + B4>) <t> + C<f> + eV<f> = 0, 



(25) 



in which the constants, A, B, C and T> are to be read as 



A = 



fa 



JO 



V = 



fa 



voR dr 



voR dr 



voR dr 



voR dr 



rt vq dR 3 
3 dr 
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a 2 2\ d-R d / tt 2 
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c rr 3 D /'d-R 

i -oR[^ 



dr, 



(26) 



respectively. The form in which equation d25t has been abstracted is that of a general Lienard system 1 StrogatJl994 ; Jordan & Smithll 999), 
All the terms of equation {25}, which carry the p arameter, e, have arisen in consequence of nonlinearity. When one sets e = 0, one readily 
regains the linear results presented earlier (Rav 2003). However, to go beyond linearity, and to appreciate the role of nonlinearity in the 
perturbation, one now has to understand the Lienard system that equation d25t has brought forth. 



5 EQUILIBRIUM AND INSTABILITY IN THE LIENARD SYSTEM 

The mathematical form of a Lienard system is like a damped nonlinear oscillator equation, going as jStrogatJl994l : ljordan & Smithlll999h 

4> + M{<j>, 4>)j> + V'(<f>) = 0, (27) 

in which, H is a nonlinear damping coefficient (the retention of the parameter, e, alongside T-i, attests to the nonlinearity), and V is the 
potential of the system (with the prime on it indicating its derivative with respect to <f>). In the present study, 



H{(t>,(t>)^A4> + B<j}, 
and 



(28) 
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V(^)=C^-+eV^-, (29) 

with the constant coefficients, A, B, C and T> having to be read from equations l |26t . 

To investigate the properties of the equilibrium points resulting from equation l !27t . it will be necessary to decompose this second-order 
differential equation into a coupled first-order dynamical system. To that end, on introducing a new variable, ip, equation j27\ can be recast 

as 



i Jordan & Smith! 19991) 



4> = ip 

ij> = -e(Att> + Bi>)1> - (Ccp + eV<f> 2 ) . (30) 

Equilibrium conditions are established with <j> = ip = 0. For the dynamical system implied by equations OOt . this will immediately lead to 
two equilibrium points on the <p-ip phase plane. Labelling the equilibrium points with a ★ superscript, one can easily see that ((/>*, ip*) = (0, 0) 
in one case, whereas in the other case, (<p*,ip*) — (—C/(eD), 0). In effect, both the equilibrium points lie on the line, ip = 0, and correspond 
to the turning points of V(<j)). Higher orders of nonlinearity may simply have the effect of proliferating equilibrium points on the line, ip — 0. 
For the present case of second-order nonlinearity, one of the equilibrium points is located at the origin of the <p-ip phase plane, while the 
location of the other will depend both on the sign and the magnitude of C /T>. 

Having identified the position of the two equilibriums points, the next task would be to understand their stability. To do so, both 
equilibrium points are to be subjected to small perturbations, following which, a linear stability analysis will have to be carried out. The 
perturbation scheme on both (p and ip is <p — <p* + 5<p and ip = ip* + Sip. Applying this scheme on equation ( 130b . and then linearising in Sep 
and Sip, will lead to the coupled linear dynamical system, 

1 {Sep) = Sip 

^ (5V) = -V"(<p*)5<t> - eH(cb*, 1>*)S1>, (31) 

in which V" (</>*) = C + 2eD<p* . Using solutions of the type, Sep ~ exp(oji) and Sip ~ exp(ajt), in equations J31l l. the eigenvalues of the 
Jacobian matrix of the dynamical system follow as 



a; = _e ? ± V e2 x" v " ( ^ ) ' (32) 

with H = H(4>* , ip*) having to be evaluated at the equilibrium points. 

Once the eigenvalues have been determined, it is now a simple task to classify the stability of an equilibrium point by putting its 
coordinates in equation d32t ■ The equilibrium point at the origin has the coordinates, (0, 0). Using these coordinates in equation d32t . one 
gets the two roots of the eigenvalues as ui = ±iVc. it can be readily seen that if C > 0, then the eigenvalu es will be purely ima ginary 



quantities, and consequently, the equilibrium point at the origin of the <p-ip plane will be a centre-type point l Jordan & Smith 19991) . And 



indeed, when the stationary inflow solution, about which the perturbation is constrained to behave lik e a standing wave, is sub-critical over the 
entire region of the spatial integration, then C > 0, because in this situation, v 2 < /3 2 Oq jRavll2003l) . Therefore, the centre-type equilibrium 



point at the origin of the phase plane indicates t hat the stan ding wave will be purely oscillatory in time, with no change in its amplitude. This 
very conclusion was made in a previous study lRavll2003h on the linearised analysis of the standing wave, and it could be arrived at equally 



correctly by setting e = (the linear condition) in equation d32t . 

The centre-type point at the origin of the phase plane has confirmed the results known already. It is the second equilibrium point that 
offers some novelties. This equilibrium point is entirely an outcome of taking nonlinearity to its lowest order (the second order) in the 
standing wave. The coordinates of this equilibrium point in the phase plane are (—C/(eD), 0), and using these coordinates in equation \32l , 
the eigenvalues become specified as 



AC , 1{AC X ' 



~ 2V ± i{zD) +C - (33) 

Noting as before, that C > 0, and that A, C and T> are all real quantities, the inescapable conclusion is that the e igenvalues, uj, are real quan- 
tities, with opposite signs. In other words, the second equilibrium point is a saddle point |jordan& Smith! 1999I) . and as such its implications 
may be far-reaching when it comes to evolving a critical solution in the axisymmetric flow through the dynamic process. 

To understand this fact, the first thing that has to be realised is that if the magnitude of the temporal part of the perturbation exceeds a 
certain critical value, i.e. if \tf>\ > \C/T>\, then the perturbation will undergo a divergence in one of its modes, i.e the stationary, sub-critical 
background solution will become unstable under the influence of the perturbation. This is how it must be in the vicinity of a saddle poin t, and 
higher orders of nonlinearity (third order onwards) will not be able to make this effect disappear (Strogatz 1994; Jordan & Smith 1999). The 
best that one may hope for is that the instability may grow in time till it reaches a sat uration level imposed by a higher order of n onlinearity, 



a feature that has a precedence in the laboratory fluid problem of the hydraulic jump JVolovildliooi iRav & Bhattacharie e 2007b). 

While all of this gives the perturbative perspective, the implications of the saddle point for the non-perturbative evolutionary dynamics 
are also noteworthy. It is evident that there can be no critical solution without gravity driving the infall process. So from a dynamic point of 
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view, gravity starts the evolution towards the critical state from an initial (and arguably nearly uniform) sub-critical state. If, however, during 
the evolution in real time, a saddle point is to be encountered, then there should be difficulties in reaching a stable and stationary critical end. 

Under linearised conditions, the perturbation on globally sub-critical flows maintains a constant amplitude. Viewed in the phase portrait, 
this feature translates into close d phase trajectories around a centre-t ype point. Now, from dynamical systems theory, centre-type points are 
known to be "borderline" cases dStroga tz 1994; J ordan & Sm ith 1999). In such situatio ns, the linearised treatment will show apparently stable 
behaviour but an instability may emerge immediately on accounting for nonlinearity i Strogatz|[l994 ; Jordan & Smith|[T999l) . This instability 
in the nonlinear analysis may be connected to the constant distribution of angular momentum in the inviscid accretion disc. To consider an 
analogous situation in an inviscid and incompressible Couette flow (which also has an axial symmetry, just like the inviscid accretion disc), 
if Rayleigh's criterion for stability is invoked, then t he stratification of an gular momentum in the flow is stable if and only if it increases 
monotonically outwards, i.e. have a positive gradient jchandrasekhaj[l981 ). If, however, the gradient of the angular momentum is negative, 
then the rotational flow will be unstable. To carry this analogy over to the inviscid accretion disc, the gradient of the constant distribution 
of angular momentum is zero. So effectively this implies a borderline case between a stable positive gradient, and an unstable negative 
gradient. Such borderline cases may show the apparently stable features of a centre-type point in a linearised dynamical sy stems analysis , 
but in these situations it is always proper to draw a definitive conclusion regarding stability only from a nonlinear analysis (Strogatz il 19941 : 
Jordan & SmithTl 999). And in the present study the nonlinear analysis does reveal an instability. A cautionary note that may be sounded here 



is tha t the Couette flow in this discussion is mathematically modelled to have an infinite extent along the vertical direction jchandrasekhaj 



1981), while the inviscid disc in this study has a finite thickness, and is modelled mathematically by taking vertically integrated quantities on 



the equatorial plane of the disc. 



6 CONCLUDING REMARKS 



The Lienard system derived in this work indicates that the number of equilibrium points will depend on the order of nonlinearity that may 
be retained in the equation of the perturbation. Additional equilibrium points, resulting from higher orders of nonlinearity, may temper the 
instability that has been found here. However, up to the second order at least, an instability in real time appears to be an undeniable fact, and it 
has been suggested that this instability, arising from the nonlinear order, could be connected to the gradient of the angular mome ntum. There 



are ot her disc models, in which angular momentum maintains a positive gradient, as for example the Keplerian accretion disc I Fran k et al 



2002). Examining the stability of such configurations, may lead to a general understanding of the connection between the stability of rotational 
flows and the distribution of their angular momenta. 

Now, real fluid flows are also influenced by viscosity. In fact, fluid flows are usually affected both by nonlinearity and viscosity, oc- 
casionally as competing effects. In models of accretion discs, viscous dissipation usually brings about stability, but in one of the proposed 
models of axisymmetric accretion, the quasi- viscous accretion disc, viscosity is see n to destabilise the flow |Bhattachariee & Ray 2007 ; 
Bhatt achariee et al. 20091) . In this model, viscosity, quantified by the a paramter of Shakura & Sunvaev 1 1973 ), is made to behave as a 
small first-order perturbative effect about a background inviscid flow. This instability, known as secular instability, is not without its prece- 
dence. ExactIy_jWsJdnd_of instability is also seen to grow in Maclaurin spheroids on the introduction of a kinematic viscosity to a first 
order jchandrasekhaj[l987 ). So when it comes to stability, the role of viscosity is not as clearly defined as it is for enabling the infall process 
by an outward transport of the angular momentum. In any case, beyond simple hydrodynamics, the stability of accretion discs is more likely 
to be governed by phenomena like radiative processes, turbulence and magnetohydrodynamics i Balbus & Hawlevl 19981) . 



As a matter of regular practice, the stability of fl uids is also studied by constraining a perturbation to behave like a travelling wave JPetterson et al 
198dlCrossll98d : lRavl2003l ; lRav & Bhattacharieel2007bl) . At ti mes, one encounters the surprising situation of a fluid flow being stable under 
one type of perturbation, but unstable under the effect of another l Cross & Hohenbergll 1993c Ray & Bhattacharje e 2007b). With nonlinearity 
lending an additional aspect, these effects are worth a close look in future studies. 
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